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We study theoretically the groundstates of two-dimensional Bose-Hubbard models which are frus- 
trated by gauge fields. Motivated by recent proposals for the implementation of optically induced 
gauge potentials, we focus on the situation in which the imposed gauge fields give rise to a pattern 
of staggered fluxes, of magnitude a and alternating in sign along one of the principal axes. For 
a = 1/2 this model is equivalent to the case of uniform flux per plaquette = 1/2, which, in the 
hard-core limit, realizes the "fully frustrated" spin-1/2 XY model. We show that the mean-field 
groundstates of this frustrated Bose-Hubbard model typically break translational symmetry. We 
introduce a general numerical technique to detect broken symmetry condensates in exact diago- 
nalization studies. Using this technique we show that, for all cases studied, the groundstate of the 
Bose-Hubbard model with staggered flux a is condensed, and we obtain quantitative determinations 
of the condensate fraction. We discuss the experimental consequences of our results. In particular, 
we explain the meaning of gauge-invariance in ultracold atom systems subject to optically induced 
gauge potentials, and show how the ability to imprint phase patterns prior to expansion can allow 
very useful additional information to be extracted from expansion images. 



I. INTRODUCTION 

One of the most striking aspects of the physics of Bose- 
Einstein condensed systems is their response to rotation. 
The rotation plays the role of a uniform magnetic field, 
which frustrates the uniform condensate, forcing it into 
a state containing quantized vortices and carrying non- 
vanishing currents^— Theory shows that at sufficiently 
high vortex density this frustration can lead to the break- 
down of Bosc-Einstcin condensation, and the formation 
of a series of strongly correlated quantum phases which 
can be viewed as bosonic analogues of the fractional 
quantum Hall states^ 

In typical magnetically trapped Bose gases^ practical 
limitations on the rotation rate (vortex density) are such 
that strongly correlated phases are expected only at a 
very low particle density where the interaction energy 
scale is very small £ As a result, it has proved difficult 
to reach this strongly correlated regime. (However, see 
Ref. [B| for interesting recent results for systems with small 
particle numbers.) 

It has been proposed that one can exploit the strong in- 
teractions that are available in systems of bosonic atoms 
confined to optical lattices^ to enhance the possibility of 
achieving these correlated phases. In this context, the 
natural model to consider is the Bose-Hubbard model 
with uniform effective magnetic flux [Eq. (fl))]. This 
"frustrated" Bose-Hubbard model can show very inter- 
esting physics, far beyond the physics of the usual Bose- 
Hubbard models Atomic systems well-described by this 
frustrated Bose-Hubbard have been studied experimen- 
tally by using rotating optical lattices*^ albeit so far 
limited to situations of large lattice constants and large 
numbers of particles per lattice site which are outside the 
strongly correlated regime. However, a series of theoreti- 
cal proposals^— indicate that it should be possible to im- 
print strong gauge fields on an optical lattice, and thereby 
realize a regime where interactions are strong, and with 



both the particle number per site, n, and vortex number 
per plaquette, n^, of order one. In this regime, theory 
shows that there are strongly correlated phases represen- 
tative of the continuum quantum Hall states limi t 11 ' 14 ! 16 
as well as related interesting strongly correlated phases 
that arc stabilized by the lattice itselfjil Other candi- 
dates are related to Mott physics i 18 ' 19 

Our confidence in the existence of strongly correlated 
phases of the frustrated Bose-Hubbard model relies on 
the results of large-scale numerical exact diagonalization 
studies ! 11 ! 14 ' 16 ? 17 However, these studies have found evi- 
dence for strongly correlated phases only in a relatively 
small region of parameter space (spanned by the particle 
density per site n, flux per plaquette n^, and interaction 
strength U/J). There are surely competing condensed 
phases, which can be viewed as vortex lattices that are 
pinned by the latticed An important question emerges 
from the point of view of these numerical approaches: 
How does one determine condensation in exact diago- 
nalization studies? In conventional condensed systems, 
one looks for the maximum eigenvalue of the single par- 
ticle density matrix of the groundstatej2i However, here 
the condensed states are (pinned) vortex lattices, and 
therefore break translational symmetry. As a result, one 
expects a degeneracy of the spectrum in the thermody- 
namic limit How does one quantify the degree of con- 
densation? 

In this paper we propose a powerful general numer- 
ical method that can be used to identify and charac- 
terize condensed groundstates which break a symme- 
try of the Hamiltonian. We use this to study several 
cases of interest in the context of optically induced gauge 
potentials^ Optically induced gauge potentials have re- 
cently been implemented experimentally without an op- 
tical lattice ^ 4 ' 25 These successes encourage a high de- 
gree of optimism that the related schemes on optical 
lattices*^ will also be successful. 

Motivated by the proposals of Jaksch and Zoller— and 
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Gcrbicr and Dalibard^ in this paper wc focus not on 
the case of a uniform magnetic field, but on the case of a 
two-dimensional square lattice with a staggered magnetic 
field, with a flux per plaquette of magnitude a that al- 
ternates in sign along one of the principal axes. This flux 
configuration involves much a simpler experimental im- 
plementation than the case of uniform flux. As described 
below, for the special case of a = 1/2 this is equivalent 
to the uniform flux. In this case, the model simulates 
a quantum version of the "fully frustrated" XY model. 
For other values of a, it represents a class of frustrated 
quantum spin models. Related but different staggered 
flux Hamiltonians can be generated by time-dependent 
lattice potentials as discussed in Ref. [2(| 

Based on numerical exact diagonalizations, we provide 
evidence showing that the groundstatc breaks transla- 
tional invariance, and is condensed for all flux densities 
and (repulsive) interaction strengths. We evaluate the 
condensate fraction and condensate wavefunctions from 
the exact diagonalization results. We describe how evi- 
dence for translational symmetry breaking can be found 
in measurements of the real-space and momentum space 
(expansion) profiles, and how these can be used to de- 
termine the condensate depletion. As part of this work, 
we explain some important general aspects of expansion 
imaging of systems involving optically induced gauge po- 
tentials. 



II. MODEL 

We shall study the properties of the two-dimensional 
Bose-Hubbard model subject to an abelian gauge poten- 
tial, as described by the Hamiltonian 
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The operator destroys (creates) a boson on the lat- 
tice site i, which we choose to form a square lattice; 
U describes the onsite repulsion (U > is assumed 
throughout); J is the nearest-neighbour tunneling energy. 
The Hamiltonian conserves the total number of bosons, 
N = Y^i 7 " l i = Ei^l^i- Throughout this work, we shall 
consider the system to be uniform, with N chosen such 
that the mean particle density per lattice site is n. The 
results of these studies can be used within the local den- 
sity approximation to model experimental systems which 
have an additional trapping potential. 

The fields (which satisfy Aij = —Aji) describe the 
imposed gauge potential. All of the physics of the sys- 
tem defined by the Hamiltonian (|T|) (energy spectrum, 
response functions, etc.) is gauge- invariant. Therefore, 
its properties depend only on the fluxes through plaque- 
ttes 
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FIG. 1: (a) The gauge- invariant flux through the 
plaquette, n^, is defined by = ^Yli,j A ij = 

2^ {A\2 + A-2'j, + A;u + An) (b) The simplest optically in- 
duced gauge potential to imprint on the square lattice^^ has 
an alternating pattern of fluxes of magnitude a, Eq. ((6]). 



where a labels the plaquette, and the sum represents the- 
directed sum of the gauge fields around that plaquette 
(the discrete version of the line integral) , as illustrated in 
Fig.QJa). Since each phase Aij is defined modulo 2-7T, the 
gauge invariant fluxes ([T]) are defined modulo 1 (i.e. are 
invariant under n a ^ — > + 1), so they can be restricted 
to the interval -1/2 < n? < 1/2. 

The gauge-invariant fluxes through the plaquettes lead 
to an intrinsic "frustration" of condensed (superfluid) 
phases on the Bose-Hubbard system. This is best under- 
stood in the case of strong interactions, U 3> J, where 
double occupancy is excluded. In this hard-core limit, the 
Bose-Hubbard model is equivalent to a spin- 1/2 quan- 
tum magnet, using the standard mapping sf = n, — i, 



§t = aj, s { = hi, with Hamiltonian (up to a constant 



shift in energy) 



h-c — 



(2) 



(The conservation of particle number becomes conser- 
vation of S z = J2i §i = & - 1/2.) This Hamiltonian 
describes a quantum spin-1/2 magnet, experiencing XY 
nearest neighbour spin exchange interactions. These ex- 
change interactions are "frustrated" by the gauge fields. 
The "frustration" can be seen by considering the natural 
mean- field limit of the spin Hamiltonian @, generalizing 
from spin-1/2 to spin-S' and taking the S — > 00 limits 
Then, the (vector of) spin operators §i can be replaced 
by the classical vector s of fixed length S. It is convenient 
to parameterize this vector as 



s = S(sin 9 cos </>, sin 9 cos (j>, cos I 



(3) 



which represents the spin by the polar and azimuthal an- 
gles 9, (f>. The Hamiltonian becomes the (classical) energy 
functional 

H m{t = -2JS 2 sin Oi sin 0j cos^ - fa + A tj ) (4) 

and it is natural define the fractional occupation of the 
lattice sites by n, = (1/2) (1 + cos£?i). 
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FIG. 2: Illustration of the frustrated coupling of the XY spin 
model around a plaquette with non-zero flux n$. For a given 
phase 0i at site 1, one can choose (f>2, 4>3, and <j)4 to maximally 
reduce the energy on three of the bonds (with a contribution 
of — JS 2 to the energy for each). However, the exchange en- 
ergy on the remaining bond is — JS 2 cos(Ai2 + A23 + A34 + 
An) = — JS 2 cos(27T7i0). For n$ 7^ 0, this bond cannot also 
be fully satisfied, indicating the magnetic frustration. Maxi- 
mal frustration occurs for n$ — 1/2. 



The restricted space of configurations with 9{ = ir/2 is 
important for groundstate configurations with (uniform) 
density n = 1/2. In this case, s z = 0, so all spins lie 
in the x-y-plane and the Hamiltonian becomes that of 
the frustrated XY model^ Here, "frustration" refers to 
the fact that, with 7^ for any plaquette, the angles 
<pi around this plaquette cannot be chosen in order to 
maximally satisfy the XY exchange couplings. This is 
illustrated in Fig. 

Maximal frustration occurs for = 1/2. This situ- 
ation is referred to as the "fully frustrated" case. The 
"fully frustrated" classical XY model has been much 
studied as an interesting frustrated classical magnet with 
a non-standard thermal phase transition^ The frus- 
trated Bosc-Hubbard model that we study here is a class 
of frustrated spin- 1/2 quantum systems that are analo- 
gous to this frustrated classical model. We shall focus on 
the nature of the groundstates of the system. 



III. STAGGERED FLUXES 

Motivated by proposals for optically induced gauge 
potentials'^ we consider a situation in which the flux is 
staggered in the x direction. Specifically, we choose the 
gauge (and geometry) described in Ref. Il5l for which 



Aij = 2irayi(xi 



Xi ) x (-iy 



(5) 



where (xi,y%) is the pair of integers that define the posi- 
tion of lattice site i (i.e. the cartesian co-ordinates po- 
sition in units of the lattice constants, a x , a y , in the x- 
and y-directions, as shown in FigQJb)). Since the lat- 
tice is square with nearest-neighbour hopping, this has 
the effect that hopping in the y direction has no gauge 
field (Aij = 0) while hopping in the x direction involves 
a phase Ajj = 2irayi(—l) Xi that alternates from one row 
to the next. The magnitude of the flux per plaquette is 



where the position x a is defined by the ^-position of the 
"bottom-left" corner of the plaquette (i.e. the minimum 
Xi for all sites i surrounding the plaquette a). The flux 
is staggered in the x-direction, as shown in Fig.QJb). 

A special situation arises for a = 1/2. Then, the case 
of alternating fluxes of = ( — l) Xa ^ is gauge- equivalent 
to that of uniform flux n% = |. The (gauge-invariant) 
properties of this case have higher translational symme- 
try than the case of a 7^ 0,1/2. Furthermore, a gauge can 
be chosen in which e lAij is real, meaning that the Hamil- 
tonian is time-reversal symmetric. Indeed the gauge ([5]) 
has this property. While gauge invariance allows the 
physics of ([T]) to be studied in any gauge, as we shall 
describe below, the expansion images of the atomic gas 
are gauge-dependent. We shall therefore be clear to spec- 
ify the gauge considered. 

Under the conditions that we study, significant insight 
into the physics of the frustrated Bose-Hubbard model 
can be obtained by studying its properties within mean- 
field theory. Indeed, one important goal of this work is to 
show how the mean-field groundstates emerge from the 
results of exact diagonalization studies. 



A. Single Particle Spectrum 

We first ignore interactions, and study the single par- 
ticle energy eigenstates. For the gauge field we consider 
(O , the unit cell can be chosen to have size 2 x 1 (in the 
x- and y-directions) , and the energy eigenstates are then 



tpi 



^(kxXi+ky-yi) 



Ipe 



Xi even 
Xi odd 



(7) 



with the momenta in the ranges —tt/2 < k x < tt/2 and 
— 7r < k y < 7r. [We express k x and k y in units of l/a x 
and l/a y respectively.] The energy eigenvalues E and 
cigenfunctions within the unit cell, (ip e ,^p ), follow from 



2J , cos(ky) cos(k x ) 

cos(k x ) cos(fc. y + 2-Kct) 



Tpo ) \lpo 

(8) 

The lowest energy state has k x = 0. For a < 
— arccos ((V5- l)/2) ~ 0.288, there is a single minimum 
at k y = ira. However, for a > 0.288 this minimum splits 
in two, and there are two degenerate minima defining sin- 
gle particle states ipf and ipf . The wavevectors of these 
states, k^' B are shown in Fig. [3l They are related by 



ky+ky = -2na. 51 Hence cos(k y 1 ) = cos(k y 1 +27ra), 
so the Hamiltonian (jHJ is the same for both states pro- 
vided the two sites of the unit cell are swapped. There- 
fore, ip A and ip B are exactly degenerate, and their wave- 
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functions related by rp^. 

The appearance of two degenerate minima in the 
single-particle spectrum leads to the question of 
whether the groundstate is a "simple" condensate, or 
"fragmented".— Following the general result,— one ex- 
pects that weak repulsive interactions will lead to a sim- 
ple condensate in which all particles condense the same 
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FIG. 3: The single-particle groundstate is non-degenerate for 
a < 0.289, but is two-fold degenerate for a > 0.289. The two 
minima are at wavevectors ky ,ky shown above (both have 
k x = 0). 



single particle state. The nature of this condensate can 
depend on the properties of the two states, and the na- 
ture of the interactions. However, if the condensate 
wavefunction contains non-zero weights of both states, 
with wavevectors ky and ky , then the condensate breaks 
translational symmetry in the y direction. 

Similar physics - of a two-component Bose gas - has 
been discussed recently in the context of optically in- 
duced gauge potentials in the continuum<2i~— For the 
case we consider, one principle difference is that there 
is an underlying lattice periodicity. Phases with broken 
translational symmetry can therefore "lock" to the lat- 
tice periodicity, leading to condensed groundstates with 
a broken discrete translational symmetry. Furthermore, 
the unit cells of the phases that we describe contain par- 
ticle currents, which can be viewed as arrays of vortices 
and antivortices. 

In the following we shall explore the nature of this sym- 
metry breaking. A principle result will be to show that 
the results of exact diagonalization studies are consistent 
with simple condensation in a symmetry broken state. 
We shall focus on two cases: 

(i) a = 1/2. Here, the model is equivalent to the case of 
uniform fluxes n a ^ = 1/2. The (gauge-invariant) system 
therefore enjoys the full translational symmetry of the 
lattice. The minima are spaced by Ak y = it, suggesting 
that the broken symmetry state will have translational 
period 2 in the y-direction. The symmetry broken state 
will have two degenerate configurations. 



0.389. This 



(ii) a = i arccos f| y |(9 — \/65) 

value is chosen such that Ak y = 2ir/3, such that the 
broken translational symmetry state has period 3 in the 
y-dircction. The symmetry broken state will have three 
degenerate configurations. 

Throughout this work, we concentrate on cases where 
the particle density n is non-integer. Thus, we do not 
discuss the properties of the Mott insulator states. While 



frustration by the gauge fields can affect the stability of 
the Mott insulating stateS ) 34 i 35 the qualitative form of 
these insulating states is unchanged. We focus on the 
nature of the superfluid states, the properties of which 
are much changed in the presence of the frustrating gauge 
fields. 



B. Gross-Pitaevskii Theory 

One can understand the effects of weak interactions, 
i.e. U <C nJ, within Gross-Pitaevskii mean- field theory. 
One assumes complete condensation, forming a many- 
body state 



N 



l* c )= 5>?al |o) 



(9) 



where ipf is the (normalized) condensate wavefunction. 
The condensate wavefunction is determined by viewing 
(0) as a variational state, and minimizing the average 
energy per particle 
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(i\r-i)X>. 
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2N 



(11) 



For large systems N 3> 1 the typical interaction energy 
(ITT1) is of the order of U n, while the kinetic energy (fTU)) 
is of order J. 

In the weak coupling limit, nil <C J, the ground- 
states of the Gross-Pitaevskii equation can be obtained 
by assuming that the condensate consists only of those 
states that minimize the single-particle kinetic energy. 
For a < a c , there is a single minimum, so the conden- 
sate will be formed from this state alone. For a > a c , 
there are two degenerate minima, and we should write 



(12) 



which (schematically) denotes a linear superposition of 
the states in these two minima. Then one chooses A and 
B to minimize the interaction energy iV'f | 4 - 



1. The case a = 1/2 

This is the case of the fully frustrated Bosc-Hubbard 
model, a = 1/2, which is gauge equivalent to uniform 
plaquette fluxes of = 1/2. Minimizing the interac- 
tion energy within states of the form (fT2"]) leads to two 
solutions which we can denote schematically by 



V2 



(13) 
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FIG. 4: Illustration of the two (degenerate) mean-field 
groundstates of the fully frustrated model, n<f, = 1/2. The ar- 
rows on each lattice site represent the phase of the condensate 
wavefunction (the amplitude is constant). The two phases are 
characterized by staggered circulating currents, in directions 
illustrated within the plaquettes. Double lines mark links 
where a 'phase' e l7r is imprinted. 



In detail, the wavefunctions are 



±.i 



2^2 - V2 



i±i(v^-i)(-i) w 
(v^ -!)(-!)* ±i 



Xi even 
Xi odd 



(14) 

where i refers to the labelling in Fig. [2] The wavefunc- 
tions are illustrated in Fig. The pattern of phases 
is equivalent to that of the ordered groundstate of the 
fully-frustrated classical XY model^. 

Since the condensed state is a superposition of states 
with different k y , it is a state with broken translational 
invariance. The difference in wavevectors is Ak y = ky — 
ky = 7T, so the new unit cell in the y-direction has size 
Ay = 2 (in units of the lattice constant a y ). 

It is useful to recall that the Hamiltonian at a = 1/2 
is time-reversal invariant, so all of its eigenvectors can be 
chosen real. The fact that the condensed wavefunctions 
(|14jl are imaginary shows that these break time-reversal 
symmetry. The two states are time-reversed partners, 

The two states are the two (time-reversed) partners of 
a so-called "staggered flux" phase, a related version of 
which has been discussed in the context of the cupratc 
superconductors^ This phase is characterised by circu- 
lating currents around the plaquettes, which are arranged 
in a staggered "checkerboard" pattern. It is straightfor- 
ward to show that the states tp± do carry these staggered 
gauge-invariant currents. We emphasize that the under- 
lying Hamiltonian is both translationally invariant and 
time-reversal invariant, so the emergence of these stag- 
gered flux states is a result of the breaking of both of 
these symmetries. 

The same condensed states have been shown to 
describe the groundstates of a related Bose-Hubbard 
models In constrast to the model that we study, in the 
model of Ref. HH, the staggered flux a is applied directly 
in a checkerboard pattern. However, for staggered flux 
of a = 1/2 the two models coincide (up to a choice of 
gauge, as described below). 



The derivation of the states (jT3|) described above was 
presented in the weak-coupling limit, nU <C J. How- 
ever, note that the state has the special property that 
the density is uniform I'i/'fl 2 = constant. Therefore, 
these condensates already minimize the interaction en- 
ergy Tlii l^i | 4 a t fixed average density. With increasing 
interaction strength nU ~ J, the condensate does not 
change. (We have confirmed this with extensive numer- 
ical simulations of the mean-field theory^) Even in the 
regime U 3> J where Gross-Pitaevskii theory becomes 
unreliable owing to the suppression of number fluctua- 
tions and a Gutzwiller mean-field theory is required, the 
condensate wavefunction remains the same. 



2. The case a = 0.389 



Choosing a = i arccos |(9 
leads to the situation in which Ak v 
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0.389 



v y — 2n/3. Then, a 
state of the form (fT2l) with non-zero A and B breaks 



translational symmetry with a period of Ay = 3 in the 
y-dircction. 

Minimizing the interaction energy over states of the 
form (fT2"]) leads to the conclusion that the optimum con- 
densate wavefunction is the superposition 



1 



</4 = -j= [V'fc.4 + e ix ip kl 



(15) 



Thus, the condensate does break translational symme- 
try, with the expected period of Ay = 3. Unlike the 
(special) case of a = 1/2, for these functions the density 
is not uniform, so spatial patterns of the density show 
this periodicity. 

One interesting feature of this mean-field result is that 
the interaction energy is independent of the phase x a P~ 
pearing in (fl~5"|) . Thus, there is an infinite set of condensed 
groundstates, which are related by the different choices of 
the phase x- In a continuum model, this infinite degen- 
eracy would correspond to the broken continuous trans- 
lational symmetry, the different choices of x denoting 
translations of the same state. Here the model is explic- 
itly defined on a lattice, so there is no continuous transla- 
tional symmetry. Since the wavefunction is defined only 
on lattice sites, the states with different values of x ar e 
not just related by simple translations. Nevertheless, for 
weak interactions nU <C J, we find that, surprisingly, 
there emerges a continuous degeneracy, beyond the ex- 
pected threefold degeneracy from translations. This gives 
an additional "Goldstonc" mode, associated with slow 
spatial variations of x as a function of position. (The 
same feature arises for values of a for which Ak y = 2n/p 
with p =/= 2.) 

Since this additional emergent degeneracy is not pro- 
tected by a symmetry of the Hamiltonian, one does not 
expect it to survive in general (e.g. to stronger interac- 
tions or to quantum fluctuations). Indeed, we find that 
by performing mean-field theory to stronger interaction 
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strengths, the continuous degeneracy is lost. Including 
energy corrections to order (nil) 2 / J we find that the 
three cases x = 0, ±2n/3 are selected as the energy min- 
ima. This gives rise to a set of three degenerate ground- 
states, with a unit cell of size Ay = 3 and all related 
by translations in the y direction. The Goldstone mode 
described above develops a gap. 



3. The case of general a 

For typical cases of a in the range a c < a < 1/2, 
there are two minima in the single particle energy spec- 
trum, but the spacing Ak y is incommensurate with the 
underlying lattice. In the weak coupling limit nil <C J, 
one expects the groundstate to remain incommensu- 
rate with the lattice, having a continuous degeneracy 
(along the lines of that described by the phase \ above) . 
However, for sufficiently strong interactions, the density 
wave will "lock" to the lattice, via an incommensurate- 
commensurate phase transition^ Additionally, for stag- 
gered flux a > a c , near the bifurcation in Fig. [3j the 
single-particle dispersion becomes very flat, and we ex- 
pect a regime of large fluctuations where potential con- 
densate groundstates will be strongly depleted. 

IV. NUMERICAL METHODS 

In order to investigate the quantum groundstate of the 
frustrated Bose-Hubbard model, we have performed large 
scale exact diagonalization studies of the model (TTJ). We 
define the system on a square lattice of N s = L x x L y 
sites, and impose periodic boundary conditions to min- 
imize finite size effects. For a total number of bosons, 
N, we therefore study a system with mean density n = 
N/(L x L y ). _ 

The possible geometries and system sizes are con- 
strained by the form of the gauge potential that is ap- 
plied. Furthermore, the plaquette fluxes determine the 
translational symmetries of the Hamiltonian, and hence 
the conserved momenta. For staggered flux ©, in the 
general case (a ^ 0,1/2) the translational symmetries 
are those of the unit cell with size 2 x 1 (in x and y direc- 
tions). Higher symmetry arises for the cases a = and 
a = 1/2, for which the flux is uniform; in the latter case, 
the (magnetic) translational symmetries of many-particle 
systems follow from Ref. [H. 

By definition, the Hamiltonian commutes with the uni- 
tary transformations that effect these translational sym- 
metries. The energy eigenstates are therefore also eigen- 
states of these translational symmetry operators (i.e. 
eigenstates of conserved lattice momentum). In any ex- 
act numerical calculation, the energy eigenvectors will 
also be eigenstates of the conserved lattice momentum. 
However, as described above, in general the mean-field 
states break the translational symmetry. Thus the mean- 
field states are not eigenstates of momentum. In order 



to make comparisons with the mean-field states, and the 
possibility of condensation, one must allow for this break- 
ing of translational symmetry. 



A. Condensate Fraction 

Often the groundstates of (repulsive) interacting 
bosons can be understood in terms of Bose-Einstein con- 
densation. Interactions between the particles lead to de- 
pletion of the condensate. If the effects of interactions 
are very strong, these may even drive a phase transition 
from the condensed phase into non-condensed phases. It 
is therefore very important to know if the groundstate 
remains condensed. (If not, then the system may be de- 
scribed by a novel, uncondensed, and possibly strongly 
correlated phase of matter.) 

The condensate fraction is quantified by the general 
definition introduced by Yangi^i From the many-particle 
groundstate l^o); one forms the single particle density 
operator 

ftj = <*o|3}a 3 -|*o>, (16) 

a Hcrmitian operator, the trace of which ^\ pa is the 
mean (total) number of particles, N. Then, one finds the 
eigenvalues of /Oy. For "simple" BECs^ the spectrum 
has one eigenvalue which is of order and which is 

therefore much larger than all others for large N (the 
thermodynamic limit). Denoting this largest eigenvalue 
Ao, the condensate density n c and condensate fraction x c 
for average density n arc defined by2i 



The eigenvector of pij corresponding to the largest eigen- 
value is the condensate wavefunction, ip®. 

While this method is applicable in the simplest of sit- 
uations, it gives misleading results in cases where the 
groundstate of the system breaks a symmetry of the 
Hamiltonian in the thermodynamic limit. In that case 
it is well known that, for a finite sized system (in which 
the ground state is an eigenstate of all symmetry op- 
erators), an analysis of the density matrix states shows 
a "fragmented" condensate in which there is more than 
one eigenvalue of order N. (See Ref. [29| for a discussion 
of fragmentation in the context of cold atomic gases.) 
The origin of this fragmentation, and its relationship to 
symmetry breaking, is well-understood for simple model 
systems with no condensate depletion, for example for 
condensation of N bosons in two orbitalsj 29 ' 40 From the 
practical point of view of exact numerical calculations, 
it is important to have a prescription for how to quan- 
tify the degree of condensation in general - i.e. in cases 
where there are many degrees of freedom and condensate 
depletion can be significant. 
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B. Condensate Fraction with Symmetry Breaking 



We propose a method to determine, on the basis of 
numerical exact diagonalization studies, the condensate 
fraction in cases where the condensed state breaks a sym- 
metry in the thermodynamic limit. Given the context of 
this paper, we focus on the case of translational symme- 
try. The method, however, is very general: no specific 
knowledge is required of the condensed state, or indeed 
of the symmetry that is broken. These emerge directly 
from the numerical calculations in an unbiased way. 

The starting point is to determine the energy spec- 
trum, in order to identify if the groundstate may have 
a broken symmetry in the thermodynamic limit. As 
ever with numerical studies, one must study the spec- 
trum with varying systems size (up to as large a sys- 
tem size as can be achieved), in order to glean informa- 
tion about the properties of the spectrum in the ther- 
modynamic limit. It is well-known that the signature of 
(translational) symmetry breaking in a finite size calcula- 
tion is the appearance of a set of quasi-degenerate energy 
levels in the spectrum, with different eigenvalues of the 
conserved quantity associated with the symmetry (i.e. 
momentum in the case of translational symmetry break- 
ing). In the context of cold atomic gases, this has been il- 
lustrated for spin-rotational invariance ) 41 ' 42 translational 
and rotational symmetry breaking ^ 22 ' 43 and parity 40 i 44 
In the thermodynamic limit, these states become degen- 
erate, and it becomes valid to superpose the states (e.g. 
as selected by an arbitrary weak symmetry breaking per- 
turbation). Any superposition of all these states forms a 
groundstate which has broken symmetry. Before super- 
position, each of the states (eigenstates of the symmetry 
operators) can be viewed as fragmented condensates^ 

Based on the results of these calculations of the energy 
spectrum, one can look to see if, in the thermodynamic 
limit, several states are tending to become degenerate. 
The emergence of this degeneracy appears when the sys- 
tem size is sufficiently large; if the degeneracy is not well 
resolved, then this suggests that the numerical calcula- 
tions are not on a sufficiently large system size to be con- 
clusive. In many cases, an emergent quasi-degeneracy 
can be very convincingly established4 i 22 ' 40 i 44 In many 
practical cases of interest where the degeneracy is not 
fully established, it can still be of value to make the hy- 
pothesis that a small number of low-energy states will be 
degenerate in the thermodynamic limit, and to test if this 
hypothesis is borne out by a high condensate fraction. 

Suppose that an analysis of the energy spectrum sug- 
gests that there are D such states, |Wq), with fi = 
1,2, . . . , D, which tend towards degeneracy in the ther- 
modynamic limit. We assume that these D states can be 
distinguished by eigenvalues of symmetry operators (e.g. 
no two have the same momenta). In order to investi- 
gate the possibility of simple BEC in a broken symmetry 



state, we propose that one forms the superposition state 

D 

I^E^o) (18) 

which depends on the D complex amplitudes c^. Then, 
for this superposition state - which is not an eigenstate 
of the (translational) symmetry - one should determine 
the single particle density matrix (fTHJ) and find the con- 
densate fraction (| 1 T[) , each of which are functions of the 
parameters c M . We define the condensate fraction of the 
broken symmetry state by in terms of the optimal choice 

X c = max[a; c (cu)] . (19) 

The corresponding optimizing coefficients define the 
associated condensate wavefunction (p~8|) . Since this is 
a broken symmetry state, in general there are D sets of 
coefficients c M which give the (same) maximum conden- 
sate fraction, and hence D such condensed states. These 
correspond to the D broken symmetry states, and are 
related by applications of the symmetry operations. 

Below, we illustrate the application of this approach 
for the cases of the Bosc-Hubbard model with staggered 
flux, at a = 1/2 (fully frustrated) where D — 2, and at 
a = 0.389 where D = 3. 



C. Unfrustrated Bose-Hubbard Model 

As a warm-up, and to test the quantitative validity of 
exact diagonalization in determining the condensate frac- 
tion, we study the Bose-Hubbard model in the absence of 
gauge fields (all plaqucttc fluxes vanish, = 0). In this 
case, for non-intcger particle density, it is known that the 
groundstate is condensed^ and the condensate fraction 
has been established by detailed numerical studies in- 
cluding quantum Monte Carlo^ as well as in spin-wave 
theory4£ For n = 1/2, the condensate fraction is x c = 1 
for U/J = 0, and falls to x c ~ 0.4 in the hard-core limit 

U/J ^ oo4£ 

We have used ED results on systems of up to L x x L y = 
5 x 6 to determine the condensate fraction. Consistent 
with the lack of symmetry breaking, in all cases the spec- 
trum shows a clear groundstate. (The broken gauge in- 
variance of this state will appear in a emergent quasi- 
degencracy at different particle numbers. We work at 
fixed particle number.) By forming the single particle 
density operator, and finding its maximal eigenvalue, we 
find that for hard-core bosons at n = 1/2 the condensate 
fraction is 0.433(6). The favourable comparison with the 
Monte Carlo result illustrates that, for this case, the sys- 
tem sizes amenable to ED arc sufficiently large to allow 
accurate quantitative determination of x c . 
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FIG. 5: Exact diagonalisation results for the condensate frac- 
tion x c for the maximally condensed state at n$ = 1/2, i.e. 
incorporating translational symmetry breaking: as a function 
of the interaction strength U (left); and as a function of par- 
ticle density n in the hard-core limit (right). In addition, the 
inset in the left frame shows the low- lying spectrum of the sys- 
tem for the case of hardcore interactions at half filling. These 
spectra clearly confirm the emergence of a twofold degener- 
ate groundstate in large systems. At N = 12, the splitting is 
barely discernable at the scale of the figure. 



D. Fully Frustrated Bose-Hubbard Model, a = 1/2 

We now turn to the case of the fully frustrated Bose- 
Hubbard model, a = 1/2, which is gauge equivalent 
to uniform plaquette fluxes = 1/2. To analyse this 
case, we study system sizes for which the translational 
symmetry^ is the largest, implying the largest possible 
Brillouin zone for the conserved (many-particle) momen- 
tum, and no degeneracy associated merely with the mag- 
netic translations. (This is the "preferred" case of d = 1 
in the terminology and notation of Ref. [38l ) 

In these cases, where no many-body degeneracy is ex- 
pected, the groundstate in the non-interacting system 
still remains degenerate due to the degeneracy of the 
single-particle groundstate. This degeneracy is split due 
to interactions in the system. We find an emerging 
two-fold quasi-degeneracy in the spectrum for sufficiently 
large system sizes even in the case of hardcore interac- 
tions, as shown in the inset of Fig. [5] (left frame). Follow- 
ing the procedure of ^IVBl we recognize this as a sign of 
possible symmetry-breaking with D = 2. We apply the 
prescription in ijlVBI to determine the maximal conden- 
sate fraction. 

In Fig. [5] we present the results of these calculations 
of the (maximal) condensate fraction for the fully frus- 
trated Bose-Hubbard model, both as a function of U/J 
for n = 1/2, and as a function of n in the hard-core limit 
J7/J — > oo. In all cases, we find that the groundstate 
appears to be fully condensed. The smallest condensate 
fraction is for n = 1/2 and U/J — > oo. Here, an ex- 
trapolation in I/TV to the thermodynamic limit yields a 
condensate fraction of about x c ~ 0.39. In the graphs 
of Fig. \5\b) there appears to be a reduction of conden- 
sate fraction at about n ~ 0.25. Wc associate this with 



the fact that, on the lattice, the Laughlin state of bosons 
would appear for full frustration = 1/2 at density 
n = 1/2 Extrapolation of the condensed fraction in 
this case yields x c ~ 0.48(2). Although the possibility 
of Laughlin correlations may act to destabilize the con- 
densed state, the groundstate at this density is a con- 
densed (supcrfiuid) phase. In all cases, the nature of the 
groundstate that we find closely matches the predictions 
of mean-field theory. First, we can check that symmetry 
breaking is of the same form. Furthermore, the conden- 
sate wavefunction that we obtain is exactly that described 
above in Eq. (flU|) . As discussed previously, these states 
have uniform density, which makes them highly robust 
to details of the interaction potential. 

The results provide the first evidence from exact diago- 
nalizations that, under all conditions, the groundstate of 
the fully frustrated Bose-Hubbard model is condensed. 
Our model includes the possibility of particle-number 
fluctuations, and thus goes beyond previous studies using 
the picture of Josephson-j unction arrays that is based on 
phase fluctuations only41 Furthermore, our results pro- 
vide a quantitative measure of the condensed fraction. 

It is interesting to compare the results to those that 
would be obtained from a Gutwillcr ansatz. In the hard- 
core limit, the Gutwillcr mean-field state has x c = 1 — n. 
Thus, at n — 1/2 the Gutzwiller theory predicts a con- 
densate fraction of 0.5. This is close to the value we ob- 
tain from exact diagonalisation results, ~ 0.39, indicating 
that Gutzwiller theory is quantitatively fairly accurate in 
this case. At n = 1/4 the Gutzwiller theory predicts a 
condensate fraction of 0.75. The exact diagonalization 
result of ~ 0.48, shows a large quantitative discrepancy. 
As described above, we attribute this to the competi- 
tion introduced by Laughlin-likc correlations which act 
to destabilise the condensate. Another interesting obser- 
vation can be made by comparing the condensate frac- 
tions for the frustrated Bose Hubbard model at = 1/2 
with the unfrustrated zero-field case. For half filling, the 
gauge-field reduces the condensate fraction seen in our 
exact diagonalizations from about x c (n$ = 0) = 0.43 to 
Xc{n<t> = 1/2) = 0.39. For n = 1/4 on the other hand, the 
reduction is more significant, with x c (n,p = 0) = 0.66(1) 
being reduced to x c {n§ = 1/2) = 0.48 in the presence of 
the field. 



E. Staggered Flux Bose-Hubbard Model, a = 0.389 

This is the staggered flux value at which Ak y ~ 2tt/3, 
so we expect a broken symmetry state with unit cell size 
Ay = 3, and thus a groundstate degeneracy of D = 3 
in the thermodynamic limit. As described above, this is 
borne out in mean-field theory, at least for sufficiently 
strong interactions. 

The numerical results are consistent with these expec- 
tations. For large interactions, a clear three-fold degen- 
eracy appears in the groundstate. Assuming D = 3 for 
all interaction strengths leads to the condensate fraction 
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FIG. 6: Exact diagonalisation results for the condensate frac- 
tion x c as a function of interaction strength U/J, for the stag- 
gered flux with a = 0.389 where the groundstate shows three- 
fold translational symmetry breaking. Results are at half fill- 
ing, n — 1/2, calculated for N = 6 particles on a system of 
dimensions L x x L y =4x3. Small symbols show optimi- 
sations over the lowest three eigenstates. For U/J < 0.1 a 
maximally condensed state requires to also include the low- 
lying Goldstone modes (large symbols; see main text). The 
condensate fraction calculated in the hardcore limit is shown 
as a dashed line. Inset: Scaling of the condensate fraction 
with system size for different U. For the largest system with 
N = 12 particles in N 3 = 24 sites, the two available lattice 
geometries yield significantly different x c . 



shown in Fig. [6] for density n = 1/2. For small U/J 
and/or density n, the results of our analysis show that 
the groundstate is condensed in the manner predicted by 
mean-field theory. For the strongest interactions (hard- 
core interactions and n = 1/2) finite-size effects remain 
significant, and it is difficult to be sure that extrapolation 
to the thermodynamic limit will leave a non-zero conden- 
sate fraction (see inset of Fig. [5]). In finite size systems, 
the condensed wavefunction obtained from the numeri- 
cal procedure is in good qualitative agreement with the 
results of mean-field theory described above. At very 
weak interactions, the presence of the Goldstone mode 
discussed in mil B 21 is also visible. For U/J < 0.1, the 
results in Fig. [6] show a discontinuous drop in x c . This 
occurs when the splitting of the three-fold groundstate 
degeneracy (due to finite size effects) is larger than the 
energy scale which mean-field theory shows is required to 
"lock" the density wave to the underlying lattice. Thus, 
this reduction in x c at small U/J is a finite-size effect. In 
this regime, we can recover a large condensate fraction, 
x c ~ 1, by including additional levels (D > 3) to account 
for the higher degree of symmetry breaking. 



V. EXPERIMENTAL CONSEQUENCES 

The phases described above break translational sym- 
metry of the underlying lattice (with 2-fold and 3-fold 
symmetry breaking in the cases discussed in detail in 



sections [ill B 1 1 and MI B 2"|) . Following the usual expecta- 
tions for symmetry breaking states, in an experiment on 
a system with a large number of atoms, one expects that 
very small perturbations (perhaps in the state prepara- 
tion) which break the perfect symmetry of the underlying 
model will cause the system to select one of the broken 
symmetry groundstatcs. (It is also possible that domains 
will form, separated by domain walls that can be long- 
lived and survive as metastablc configurations.) 

In general, one expects that this translational sym- 
metry breaking will appear in the real-space images of 
the system (i.e. if in situ imaging on the scale of the 
lattice constant is possible). This is the case for the 3- 
fold symmetry breaking described in mil B 21 There, the 
real-space image of the particle density will have spa- 
tial structure with a unit cell that has size 3x2, and 
is therefore 3 times larger than the unit cell of the un- 
derlying microscopic model. The three broken symmetry 
phases can be distinguished by the three possible posi- 
tions of this unit cell. However, for = 1/2, where 
there is a 2-fold degeneracy of the groundstate, the two 
broken symmetry phases cannot be distinguished by the 
real-space image. Here, the groundstate is the staggered 
flux phase. This has both broken translational symme- 
try and broken time reversal symmetry. However, it is 
invariant under the combined action of translation and 
time reversal. Since density is time-reversal invariant, the 
state has uniform density. Thus, one cannot distinguish 
the symmetry breaking in real-space images. (In non- 
equilibrium situations where there can be domain walls 
separating different phases, there may appear density in- 
homogeneities associated with the domain walls.) As we 
now discuss, there are ways to detect these symmetry 
broken states in the expansion images. 



A. Expansion Images 

The nature of the groundstate can be probed by ex- 
pansion imaging. We assume that all fields are released 
rapidly (compared to the subband splitting of the lat- 
tice), and that the particles (of mass M) expand bal- 
listically - without any potentials, or gauge fields and 
neglecting further interactions - according to the free- 
particle Hamiltonian, denoted Hf vcc . Then the expansion 
image after a time t is given by-^ 



n(x) = (M/htf \w{k)\ 2 G{k) 



(20) 



where k = Mx/ht, w(k) is the Fourier transform of the 
Wannier state of the lowest Bloch band, and 



G(k) 



1 



N s ^ 



\" e ifc-(r 4 



a a,- 



(21) 



is the Fourier transform of the single particle density 
matrix. For states that are well-described as conden- 
sates (i.e. with small condensate depletion), the density 



matrix is 



id] in 



(ipi)* V*? i where tpf is the condensate 



10 



wavcfunction. Thus, the expansion image directly pro- 
vides (the Fourier transform of) this condensate wave- 
function. As described in detail below in §VB1 depletion 
leads to a reduction of the amplitudes of the "conden- 
sate" peaks in the expansion image and to the appear- 
ance of an additional incoherent background. 

Compared to the usual Bose-Hubbard model, in the 
situation that we consider here - of an optically induced 
gauge potential on the lattice - there are three important 
new considerations concerning expansion images. 

• The first aspect relates to gauge invariance. Con- 
sider making a change of the vector potential in the Bosc- 
Hubbard Hamiltonian from to 



A' iA = At 



Si Sj 



(22) 



where Si is any set of real numbers. This "gauge trans- 
formation" leaves the fluxes ([TJ unchanged, which are 
therefore said to be "gauge-invariant". The new Hamil- 
tonian (with A\j in place of Aij ) can be brought back to 
its original form by introducing the operators 



ate 



-iSi 



(23) 



In this way, the Bose-Hubbard Hamiltonian adopts the 
same form as at the start (again with gauge fields Aij), 
but now with a' replacing a. Any property that is insen- 
sitive to the distinction between a' and a in (|23j) remains 
the same, and is therefore gauge invariant. Such quanti- 
ties include the energy spectrum, density response func- 
tions, etc.; indeed any observable of the closed system 
described by ([I]) is gauge-invariant. All of these gauge- 
invariant properties depend only on the gauge-invariant 
fluxes H|). 

An important point is that the expansion image is not 
gauge- invariant. Under the transformation (|23| the sin- 
gle particle density operator becomes 



(a?a'j) 



(24) 



The gauge transformation affects the Fourier transform 
of the density operator ([21]) . and therefore the expan- 
sion image (|20|). There is no inconsistency with general 
principles of gauge invariance. As described above, prior 
to expansion, all physical properties of the Bose-Hubbard 
Hamiltonian |T]) are completely unchanged. The essential 
point is that the expansion image involves the evolution 
of the system under the free-space Hamiltonian, Hf rC c, 
and this Hamiltonian is unchanged (i.e. the gauge trans- 
formation was not applied to it). Indeed, if for the expan- 
sion imaging all optical dressing is switched off, then this 
Hamiltonian is always in a fixed gauge with vanishing 
vector potential. Provided the expansion image is taken 
under the evolution of a Hamiltonian Hf lcc with vanish- 
ing gauge potential, the expansion image measures the 
canonical momentum distribution of the particles, and 
therefore depends on the gauge used in the Bose-Hubbard 
Hamiltonian H before expansion. 

• The second consideration relates to the fact that, 
owing to the optical dressing, the atoms are in more 



than one internal state. In particular, for the schemes 
of Refs. [HI and [^, the atoms on alternating sites along 
the x axis are, in turn, in the ground state, |<jr), and an 
excited state, |e), of the atom. Therefore, upon release 
of the cloud, one has the possibility to study expansion 
images of several different types. The results described 
above (|20I21[) apply only if the image at time t is formed 
in such a way that the measurement does not distinguish 
between the different internal states of the atom. How- 
ever, since the electronic states, \g) and |e), are very dif- 
ferent, it is also possible to perform measurements which 
are state-specific: one can form the expansion image of 
the \g) atoms, or of the |e) atoms. In these cases, one 
should replace G(k) in ([2T?]) by 



G 9/^) = ^- elk<r ' 



(25) 



where the change is that the sums should be over those 
sites i,j on which the atoms are of type g or e. 

• Finally, the third consideration - although not re- 
stricted to systems of this type - arises naturally from 
the ability to address sites of type g and e separately by 
spectroscopic methods. This allows for much more inter- 
esting and useful possibilities in the expansion imaging, 
including the possibility to imprint phase patterns prior 
to expansion. Specifically, immediately prior to expan- 
sion, one can choose to drive (coherent) transfer of the 
atoms from one state to an other. For example, using 
coherent laser fields of the same type as used to pro- 
vide the laser-assisted tunneling^^ one can choose to 
transfer all atoms (initially of both e and g type) into a 
given "target" state (a superposition of e and g) while 
adding a spatially dependent phase, S^ . Alternatively, 
the phases can be imprinted by site-selective potentials 
Vi applied for a short time t, leading to S° xp = Vit/h. In 
the language of the earlier discussion, immediately prior 
to expansion one has effectively applied a gauge transfor- 
mation to the initial wavefunction. The expansion image 
follows from (|2"D|) but now with 



(a^'j) 



AST 



(26) 



replacing the density operator in (f2"Tj). (Since we chose 
the same target state for all atoms, they are indistin- 
guishable in the final image so all sites contibute.) This 
additional freedom is not restricted to dressed atomic sys- 
tems of the type we have described. Indeed, the applica- 
tion of a spatially- varying potential V(r) over a time t to 
a one-component condensate will cause the local phase 
to wind by V(r)t/h. If the time t is short compared to 
microscopic timescales, this can be viewed as an instan- 
taneous phase-imprinting prior to expansion. Techniques 
of this kind could be used to tune out the "shearing" in 
the expansion images of Ref. [25|. 

The additional freedom to imprint phases greatly en- 
hances the information that can be extracted from ex- 
pansion images. As we illustrate below, it allows impor- 
tant information regarding the nature of the phases to 
be extracted. 
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FIG. 7: The expansion image for the staggered flux phase, as determined by the structure factor as a function of (k x , k y ). (a) 
For the gauge field considered (e = 0) the two broken symmetry states have the same expansion image, (b, c) With a phase 
imprinting before expansion, denoted by e, the two states have different expansion images. 



1. Expansion Images for Full Frustration, a = 1/2. 

As explained above, the groundstate at n§ = 1/2 is the 
"staggered flux" phase, which is two-fold degenerate. In 
the gauge ([5]) we have been considering, the condensate 
wavcfunctions are given by Eq. (|14p . and therefore have 
a unit cell of size 2 x 2. A straightforward calculation of 
G ± (k) shows that the two broken symmetry states have 
the same Fourier transform, G + (k) = G~(k), illustrated 
in Fig. [7] (a) . Therefore, the expansion image cannot dis- 
criminate between whether the groundstate is in state "0+ 
or in state -01. On the other hand, the expansion image 
can distinguish these states from the expansion images 
of condensates formed from all particles condensed in ei- 
ther one or the other of the two single-particle states. 
These individual single-particle states have the full trans- 
lational symmetry of the underlying system - namely un- 
der x — > x + 2a x and y — > y + a y - so the Fourier compo- 
nents of a condensate formed from either one has peaks 
spaced by the reciprocal lattice vectors K x = (ir/a x ,0) 
and K y = (Q,2n/a y ). On the other hand, the "stag- 
gered flux phase" has broken translational invariance in 
the y direction, being invariant only under the transla- 
tions y — > y + 2a v , so the Fourier components are spaced 
by the (-Tr/a^O) and (0, ir/a y ). The appearance of this 
smaller periodicity in the k y direction is indicative of the 
broken spatial periodicity. 

For the gauge used in this work, the expansion image 
does not distinguish between the two different s tag gered 
flux states. However, in the gauge used in Ref. [26| these 
two states have very different expansion images. Indeed, 
in the gauge of Ref. [26| one state has a condensate wave- 
function with uniform density and phase (and therefore 
enjoys the full symmetry of the lattice), while the other 
has a phase pattern with a 2 x 2 unit cell; these give rise 
to very different Fourier transforms and therefore expan- 
sion images. 
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FIG. 8: A pattern of phases to imprint before expansion. For 
e / 0, this causes the expansion images of the two staggered 
flux states to differ. 



As described above, the ability to locally address sites 
of the lattice leads to the possibility to apply a spatial 
phase pattern e ; to the system prior to expansion. This 
can be used to discriminate between the two staggered 
flux states. Specifically, we consider the case in which 
potentials or optical dressing is used to imprint the phase 
pattern 



Si ~ e x mod(xi, 2)mod(j/j, 2) 



(27) 



Thus, for the atoms on sites with Xi = odd (i.e. which 
are all of "excited" or all "ground" states in the pro- 
posal of Ref. [TBh . the phase of every other one along 
the y-direction is advanced by e, as illustrated in Fig. [8] 
This can be achieved, for example, by applying a state- 
selective potential with period Ay = 2. 

A calculation of the resulting expansion images shows 
that this phase pattern causes the two staggered flux 
phases to have different expansion images, G+(fc,e) ^ 
G-(k,e). Owing to the relation — (tp^_)*, it is 
straightforward to show that G_(fc, e) = G+(fc, — e). The 
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FIG. 9: The expansion image of the staggered flux state is 
very sensitive to e / 0, here shown for Even a phase 

of e = 7r/10 allows for notable change. Since G_(fc,e) = 
G+(k, — e), this change allows for a discrimination between 
and 



results arc illustrated for e = ir/2 in Fig. E£b) and (c). 
There is a clear distinction between the expanstion im- 
ages of the two states: e.g. the signal at (k x ,k y ) = (1, 1) 
is absent for tp_, while being strong for ip+. It is not, 
however, necessary to impose large phase difference Si 
to obtain large effects. The change in the spectrum in- 
creases linearly for non-zero e. As shown in Fig. Oeven 
very small changes of phase e can give rise to notable 
changes in the expansion images. 



B. Measurement of the Condensate Fraction 

The condensate fraction can be measured experimen- 
tally by analysing the expansion images of the lattice 
gas. For a perfect condensate only few a coherent peaks 
are visible within the Brillouin zone. Condensate deple- 
tion from strong interactions in the atomic gas results in 
the appearance of an additional background, i.e. some 
of the density of particles is spread out in the Brillouin 
zone. Generally, we may write the density matrix of the 
depleted condensate as 



Pij = n c (V>i)"V| + 8pij, 



(28) 



where n c is the condensate density and rpf is its wave- 
function; this equation thus defines Spij. Similarly, the 
amplitudes of the expansion image can be written as 



G(k) = n c G c (k) + AG(k), 



(29) 



where the coherent part G c (k) derives from the non- 
interacting condensate 



G c (k) 



i,3 



k(r 



(30) 



This defines AG(fc). Experimentally, one can only mea- 
sure G(k). In order to extract n c = nx c , one needs to 
make some assumption about AG(k). Numerically, we 
find that this background of the expansion image has 
some internal structure (and this data could be used to 
build a more accurate model), but to a first approxima- 
tion we may assume that it is homogeneous. This trans- 
lates into making the simple assumption that AG(k) is 
independent of k. One then obtains 



G(k) = n c G c (k) + (n-n c ), 



(31) 



which satisfies the proper normalisation of the Fourier 

amplitudes G ( k ) = J2i Pa = N - 

Let us test the accuracy of this assumption for the 
example of a half filled lattice at a = 1/2 that was dis- 
cussed in W A II To visualize the effect of condensate 
depletion, Fig. [10] displays the evolution of the expansion 
image for the case already displayed in Fig. |7£a) for a 
weakly interacting gas. Unlike in the single particle pic- 
ture, calculations of the actual many-body wavefunction 
for the interacting system are limited to finite size. Data 
in Fig. [TU] were obtained for a lattice of size 4 x 4, so 
the 'incoherent' background occurs at peaks spaced by 
Ak x = 2tt/L x = 7r / '(2a x ) and Ak y = 2ir/Ly = ir/(2a y ). 
The main features of the expansion image remain those 
of the pure condensate even with strong interactions, ex- 
cept for the partial suppression of the coherent peaks that 
proceeds according to (|29|) with AG(k) given by small 
amplitudes. The Brillouin zone can thus be partitioned 
into areas with peaks Ap and the remaining background 
area Abg- The magnitude of this corrective term AG(k) 
is shown in Fig. [TJb). Its spatial dependency is weak, 
excepting the notably smaller correction within .Ap as 
compared to the overall background Abg- Therefore, 
the average amplitude in the background 

«-n^(I(i BG ))=^ ]T G(k) (32) 

is a proxy for the level of condensate depiction (n — n c ). 
Applied to the case with hardcore interactions in Fig. [TU1 
we deduce a condensate fraction of n c = 0.214 from the 
intensity of the background, which should be compared 
to the exact value of n c = 0.256. This estimate is rather 
crude, as the background also contributes some signal 
within .Ap. A more accurate value is obtained allowing 
for such a contribution (SI(Ap)) = k(/(^4bg))> where k is 
a 'coherence' factor for the addition between the coherent 
condensate wavefunction and the incoherent background. 
The corrected estimate becomes 



-4 B g + kA p 



A 



BG 



A* 



(I(Abg)), 



(33) 



which reduces to (|32p in the limit of sharp coherence 
peaks _4p — > 0. For the data in Fig. [TUl we find exact 
values of k in the range of 0.4 to 0.5. Assuming n = 0.45, 
the exact condensate depletion is reproduced to within 
1% accuracy. 



13 



a) 



7t/a y _ a 



Jt/2a 




« 




-7t/2a 



-Ji/2a 



b) 



O u=0.01, x c =0.999 
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FIG. 10: a) Expansion images G(k) of interacting Bose- 
Einstein condensates, showing the effect of condensate de- 
pletion. The data show the evolution of the condensate at 
ri = = 1/2 on a square lattice of size L x — L y = 4 for 
weak to hardcore interactions. The depletion of the coherent 
peaks is given approximately by the condensate fraction, b) 
The bottom panel shows the magnitude of the corrective term 
AG(k) in Eq. ((29} . This data shows some structure, in par- 
ticular, the background contribution is smaller at the fc-points 
where coherent peaks were present then it is elsewhere. 



VI. SUMMARY 



rise to a density wave pattern. For a = 1/2 the mean- 
field groundstate breaks both translational symmetry 
and time-reversal symmetry, forming the staggered flux 
phase which has uniform density but circulating gauge- 
invariant currents. 

We have introduced a general numerical technique to 
detect broken symmetry condensates in exact diagonal- 
ization studies. Using this technique we have shown that, 
for all cases studied, the Bose-Hubbard model with stag- 
gered flux a is condensed. We have obtained quantitative 
determinations of the condensate fractions. In particu- 
lar, our results establish that the fully frustrated quan- 
tum XY model is condensed at zero temperature, with a 
condensate fraction of x c ~ 0.4. 

The low-temperature condensed phases that appear 
in this system are of significant interest in connec- 
tion with their thermal phase transitions into the high- 
temperature normal phase. The groundstate of the fully 
frustrated system breaks both U(l) and Zi symmetries 
(owing to both Bosc-Einstein condensation, and the com- 
bination of translational symmetry breaking and time- 
reversal symmetry breaking). The transition to the high- 
temperature phase is interesting, combining the physics 
of the Kosterlitz-Thouless transition [?7(1)] with an Ising 
transition (Z2), and its properties have stimulated much 
theoretical debate i 48 ' 49 The cold atom system described 
here will allow this transition to be studied also in a 
highly quantum regime, with of order one particle per 
lattice site, that is inaccessible in frustrated Josephson 
junction arrays4£ Our results show that the staggered 
flux model leads also to other cases that break U(l) 
and Z3 (or higher Z p symmetry depending on the flux 
a). There is also the possibility to study commensu- 
rate/incommensurate transitions driven by locking of the 
density wave order to the underlying lattice. 

We discussed in detail the experimental consequences 
of our results. In our discussion, we have explained the 
meaning of gauge-invariancc in ultracold atom systems 
subject to optically induced gauge potentials: Expansion 
images are gauge-dependent, and (provided gauge fields 
are absent during expansion) measure the canonical mo- 
mentum distribution. Furthermore, we have shown how 
the ability to imprint phase patterns prior to expansion 
(analogous to an instantaneous change of gauge) can al- 
low very useful additional information to be extracted 
from expansion images. 



We have studied the groundstates of two-dimensional 
frustrated Bose-Hubbard models. We have focused on 
the situation in which the imposed gauge fields give rise 
to a pattern of staggered fluxes, of magnitude a and of 
alternating sign along one of the principal axes. For 
a = 1/2 this model is equivalent to the case of uni- 
form flux 7i0 = 1/2, which is the "fully frustrated" XY 
model with time-reversal symmetry. We have shown 
that, for a c < a < 1/2, with a c « 0.389, the mean- 
field groundstate breaks translational invariance, giving 
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